A general approach to few-cycle intense laser interactions with complex atoms 
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A general ab-initio and non-perturbative method to solve the time-dependent Schrodinger equa- 
tion (TDSE) for the interaction of a strong attosecond laser pulse with a general atom, i.e., beyond 
the models of quasi-one-electron or quasi-two-electron targets, is described. The field-free Hamil- 
tonian and the dipole matrices are generated using a flexible B-spline it-matrix method. This nu- 
merical implementation enables us to construct term-dependent, non-orthogonal sets of one-electron 
orbitals for the bound and continuum electrons. The solution of the TDSE is propagated in time 
using the Arnoldi-Lanczos method, which does not require the diagonalization of any large matri- 
ces. The method is illustrated by an application to the multi-photon excitation and ionization of 
Ne atoms. Good agreement with 7?-matrix Floquet calculations for the generalized cross sections 
for two-photon ionization is achieved. 

PACS numbers: 32.80.Fb,32.80.Rm,42.65.Re 



I. INTRODUCTION 

The ongoing development of ultra-short and ultra- 
intense light sources based on high-harmonic generation 
and free-electron lasers is providing new ways to gen- 
erate optical pulses capable of probing dynamical pro- 
cesses that occur on attosecond time scales These 
attosecond pulses are providing a window to study the 
details of electron interactions in atoms and molecules in 
the same way that femtosecond pulses revolutionized the 
study of chemical processes. Single attosecond pulses or 
pulse trains open up new avenues for time-domain studies 
of multi-electron dynamics in atoms, molecules, plasmas, 
and solids on their natural, quantum mechanical time 
scale and at distances shorter than molecular and even 
atomic dimensions. These capabilities promise a revolu- 
tion in our microscopic knowledge and understanding of 
matter Q. A major role for theory in attosecond science 
is to elucidate novel ways to investigate and to control 
electronic and other processes in matter on such ultra- 
short time scales. 

The ingredients of an appropriate theoretical and com- 
putational formulation require an accurate and efficient 
generation of the Hamiltonian and electron— field inter- 
action matrix elements, as well as an optimal approach 
to propagate the time-dependent Schrodinger equation 
(TDSE). Many theoretical papers have been devoted 
to the propagation of the TDSE including laser pulses. 
The earliest calculations employed finite-difference meth- 
ods [|[ to discretize the spatial coordinates. As shown 
in a recent review by Pindzola et al. this method 
is still being used with great success today. Other for- 
mulations employ finite-element [E], discrete- variable, or 
finite-element discrete- variable representation (FEDVR) 
[E> 0) [1] approaches to discretize the coordinates and 
thereby take advantage of the higher accuracy afforded 
by these methods. Time propagation of the wavefunc- 



tion may also be accomplished by a variety of techniques. 
These include simple approaches such as the leapfrog 
or Runge-Kutta [E] method to more sophisticated split- 
operator or Krylov space iterations [III E2 ■ A se- 
lected set of references is given in the bibliography. The 
relevant physical information is extracted from the TDSE 
by projecting the wavefunction onto appropriate long- 
range solutions after the laser interaction has vanished. 
The details of the process depend on what parameters 
are desired; total ionic yields are relatively simple to ex- 
tract while differential or doubly differential quantities 
necessitate more work Q. 

In this paper we consider a new approach to model 
the interaction of an atomic system with a strong 
laser pulse. We combine a highly flexible i?-matrix 
method [TU, EH, EE], including non-orthogonal sets of 
atomic orbitals to describe the initial bound state as well 
as the ejected-electron— residual-ion interaction, with the 
Arnoldi-Lanczos iterative propagation scheme. In con- 
trast to many other methods currently being used for 
such problems [IE, E3, EH j the present implementation 
is not restricted to (quasi-)one- or (quasi-) two-electron 
targets. It can be applied to complex atoms, such as in- 
ert gases other than helium and even open-shell systems 
with non-vanishing spin and orbital angular momenta. 
We illustrate the method with results for multi-photon 
excitation and ionization of neon by a linearly polarized 
laser pulse. 



II. NUMERICAL METHOD 

A. The B-Spline R-Matrix Method 

Unless specified otherwise, atomic units are used 
throughout this manuscript. The TDSE for the N- 
electron wavefunction ^(ri, . . . , rjy; i) of the present 
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problem is given by 

9 T/ 

[-ff (ri, ...,r N ) + V(rx, ...,r N ;t)]^f(n, ...,r N ;t), 

(1) 

where Ho(vi, rjy) is the field- free Hamiltonian con- 
taining the sum of the kinetic energy of the N elec- 
trons, their potential energy in the field of the nucleus 
(we assume an infinite nuclear mass), and their mutual 
Coulomb repulsion, while V{r\, Vn; t) represents the 
interaction of the electrons with the electromagnetic field. 
We expand the wavefunction as 

*(ri,...,rAr;t)=£C,(t)$ g (ri,...,r JV ). (2) 

Here $> q (ri, ...,rjv) are a known set of iV-electron states 
formed from appropriately symmetrized products of 
atomic orbitals. Optimization procedures tailored to the 
individual neutral, ionic, and continuum orbitals may be 
employed since the atomic one-electron orbitals are not 
forced to be orthogonal. The radial parts of the atomic 
orbitals are themselves expanded in S-splines. In the 
practical implementation of the B-spline i?-matrix (BSR) 
method [l5[ , factors that depend on angular and spin mo- 
menta are separated from the radial degrees of freedom. 
This enables the production of a "formula tape" since 
many Hamiltonian matrix elements share common fea- 
tures. Given a set of atomic orbitals, it is possible to 
realize a great economy in the construction of the ac- 
tual Hamiltonian matrix using this symbolic tape even 
for non-orthogonal basis sets, since ultimately every ma- 
trix element is a linear combination of one-electron and 
two-electron radial integrals multiplied by overlaps and 
angular factors. 

A significant advantage of the BSR method in the cal- 
culation of both bound and continuum states is the pos- 
sibility of employing non-orthogonal sets of atomic or- 
bitals for different target states, thereby omitting the 
need for pseudo-orbitals to account for the strong term- 
dependence that exists in many complex targets, with the 
noble gases being a prime example. Furthermore, we do 
not force the partial wave describing a continuum elec- 
tron with orbital angular momentum t to be orthogonal 
to all bound orbitals with the same I. While the method 
gives great flexibility in the target description, allowing 
for accurate representations with relatively small configu- 
ration interaction expansions, and also simplifies the gen- 
eral form of the close-coupling expansion used to generate 
the bound and excited states of the atomic system, the 
price to pay is the representation of the field-free Hamil- 
tonian and the dipole matrices in a non-orthogonal basis. 
If desired, the non-orthogonality of the primitive B-spline 
basis could easily be removed by replacing the splines 
by another complete but orthogonal basis, e.g., a finite- 
element discrete- variable representation [8j. However, if 
one wants the flexibility associated with a non-orthogonal 



set of physical orbitals expanded in any primitive basis, 
it is necessary to deal with the non-orthogonality issues 
directly. 

The interaction of the atomic electrons with the time- 
dependent electric potential, in the length form of the 
electric dipole approximation, is given by 

N 

V(r 1 ,...,r N ;t) = J2 E ( t )- r i ( 3 ) 

»=i 

where E(t) is the electric field. This form has been used 
for the calculations in this paper. For simplicity of the 
notation, we have omitted the spin-coordinates of the 
electrons. Since the initial bound state is a singlet state 
in our case, only singlet states will have to be coupled in 
the subsequent partial-wave expansion. 

The tasks at hand are now i) the preparation of the 
initial state, ii) the time propagation of the C q (t), and 
iii) the extraction of physically relevant information from 
the final state after the time propagation. As men- 
tioned above, the present appro ach employs the BSR 
method described in refs. [Rg, US, UM to compute all the 
time-independent matrix elements needed for the prob- 
lem. These parts require a representation of the field- 
free Hamiltonian matrices for the partial-wave symme- 
tries 1 S e , 1 P°, 1 P e , 1 D e , 1 D°, . . ., as well as the dipole 
matrices that couple any given value of the total orbital 
angular momentum L with a given parity to the sym- 
metries with L and L ± 1 of the opposite parity. All of 
these matrices can be readily generated with the BSR 
method, which may also be used to represent the initial 
bound state. Since the time dependence of the Hamilto- 
nian appears as a simple multiplicative factor, this only 
needs to be done once at the beginning of the calcu- 
lation. When the expansion in |(3J) is inserted into the 
Schrodinger equation, we obtain, 

8 N 
z5-C= [H {ri,...,r N )+J2E{t)-ri]C, (4) 

i— 1 

where S is the overlap matrix of the basis functions. 
Since we are initially interested in excitation and single 
ionization of the target atom by the laser pulse, the sym- 
metries of the field-free Hamiltonian must also contain a 
sufficient number of singly excited bound states as well 
as the continuum states representing electron scattering 
from the residual ion. As a method developed to treat 
exactly such problems, the BSR approach is particularly 
suitable to represent these states. 

B. Time Propagation 

Time propagation of the initial wavefunction may 
be accomplished using a number of approaches. Ex- 
plicit, norm-conserving approaches, which rely on simple 
matrix-vector multiplication, are generally preferred to 
implicit methods, which require the solution of a set of 
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linear equations. Of the former methods, we found the 
Arnoldi-Lanczos approach [ll], [l9[ to be quite effective, 
provided one is able to deal with the (often) poorly condi- 
tioned matrices generated by non-orthogonal basis sets. 
A general discussion and error analysis of the Arnoldi- 
Lanczos method can be found in the work of Saad [20l |. 
Here we only sketch the basic ideas relevant to our solu- 
tion of the TDSE. 

A straightforward approach is to transform the non- 
orthogonal many-electron basis to an orthogonal basis 
using the Lowdin transformation to generate new field- 
free Hamiltonian and dipole matrix blocks through 

H' = S-^ 2 H Q S-^ 2 , (5) 

D' = S-^DS- 1 ' 2 . (6) 

We thus have 

i^C = H'(t)C, (7) 

where H'{t) = H' + E(t)D'. Since H , £), and S are 
all time-independent, this only requires the diagonaliza- 
tion of the overlap matrix once, and a few matrix- vector 
multiplications at every time step. 

The essential idea of the Arnoldi-Lanczos method is 
to construct a reduced Krylov space of dimension m, at 
time t + At, 

K. m (H', v) = span{i;, H'v, H ,2 v, H'^'^v}, (8) 

where the initial vector v is the previously computed so- 
lution at time t. These vectors, which are generated by 
repeatedly applying the Hamiltonian H'(t) on the vec- 
tor v, are not used directly, but orthonormalized using 
the Lanczos recursion, 

Pn+lVn+l = {H' - a n )v n - P n v n -x (9) 

to transform the Hamiltonian matrix to tridiagonal 
form [20( 1 as long as the original matrix is Hermitian. 
The elements, a n and (3 n , of the tridiagonal matrix, may 
be computed (see below for a slightly more general case) 
during the recursion process using simple scalar prod- 
ucts. The resultant tridiagonal matrix is then diagonal- 
ized using standard algorithms. The result of the above 
procedure is an N x m matrix Q, which transforms the 
matrices from H' with rank N to h with rank m. Finally, 
the time evolution from t to t + At is achieved through 

C'(t + At) = Qe-* hAt Q^C'(t). (10) 

At each step m of the process, a convergence test is per- 
formed and once the propagated solution from two suc- 
cessive time steps has fallen below a predetermined crite- 
rion, the recursion is terminated. As long as the rank m 
of the process is substantially smaller than the original 
matrix size N, the process can be very effective. Finally, 
we note that the Arnoldi-Lanczos algorithm outlined 
above conserves the norm, i.e., \C'(t + At)\ 2 — \C'(t)\ 2 . 



One of the appealing features in the Arnoldi-Lanczos 
procedure is the fact that only matrix-vector multiplica- 
tions and scalar products are required. This allows us 
to take advantage of specific algorithms if the matrix is 
sparse. It is also worthwhile to make a few remarks re- 
garding the size of the Krylov space and the time step to 
obtain stable and accurate solutions. Since we need to 
generate new vectors in the Krylov space repeatedly and 
hence want to keep the size of that space manageable in 
practical calculations, we can only take relatively small 
time steps. For most calculations presented in this paper, 
1000 time steps per optical cycle were used to propagate 
the system. Not surprisingly, as already noted by Park 
and Light 11], numerical experiments showed that en- 
larging the Krylov space allows for larger time steps to 
be taken. This relationship was used to optimize the 
efficiency of our algorithm. 

An alternative, theoretically equivalent approach gen- 
eralizes the Lanczos process to a non-orthogonal basis. 
It thus allows for directly solving Eq. (fj| without trans- 
forming it to an orthogonal basis. In this case, we use 
the recursion 

f3 n+1 Sv n+1 = (H - a„S)v n - l3„Sv n -i = q n , (11) 

where the v n are the Lanczos vectors. These vectors are 
required to satisfy the condition, 

vlSv m = S ni7n . (12) 

This so-called S'-orthogonalization is possible since S is 
positive definite. The calculation proceeds along the fol- 
lowing steps. After computing 

a n = vlHv n , (13) 

the q n may be generated through matrix-vector multi- 
plication and previously obtained coefficients. The next 
step computes 

I3 n+1 = \jqlS- l q n . (14) 

In practice, no matrix inversions are performed. Instead, 
the ^-matrix is decomposed using the Cholesky decom- 
position for positive definite matrices, 

S = L^L, (15) 

at the beginning of the calculation. This yields 

n+ i = \JrXT n , (16) 

with 

T„ = (L-y qn , (17) 

and only requires the solution of a triangular set of linear 
equations once the Cholesky decomposition is performed. 
To complete the calculation it is necessary to solve a sec- 
ond set of triangular equations, namely 

v n+1 = L- l T n /p n+1 . (18) 
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From a numerical point of view, the Cholesky decomposi- 
tion is somewhat cheaper than the diagonalization of the 
overlap matrix, while the solution of the triangular sets 
of linear systems at each iteration is comparable in cost 
to the matrix-vector multiplication. For the case of an 
orthonormal basis, S — J, the algorithms are identical. 
More importantly for future work, we note that there are 
other possibilities, which entirely avoid the need for ei- 
ther inverting or decomposing any large matrices during 
the calculation (2l| . 

For the present work, we implemented the Arnoldi- 
Lanczos method with a fixed size of the Krylov space. As 
expected, results obtained with either of the above meth- 
ods were identical, and the matrix sizes we had to deal 
with (ranks of less than 500 for each individual block of 
the field-free Hamiltonian and the dipole matrix) were so 
small that we could actually check the results by perform- 
ing an exact diagonalization. A well-known alternative 
regarding the size of the Krylov space requires checking 
the convergence at each step and, if necessary, augment- 
ing the size of the space. By comparing results obtained 
with different sizes of the Krylov space, we ensure nu- 
merical convergence of the final results with the size of 
that space. 

In addition, we check the dependence of the final re- 
sults on the number of coupled symmetries. It is well 
known [22| that the number of L-values to couple in- 
creases strongly with decreasing laser frequency. Em- 
ploying the velocity gauge to express the dipole opera- 
tor is expected to reduce this problem [23j. However, 
the velocity gauge is problematic at short distances, and 
the pro blems increase with the nuclear charge of the tar- 
get [24j . Hence, we plan to explore switching between 
gauges at some distance in future work. Finally, to en- 
sure that box effects do not disguise the actual physics, 
we also use a masking function to avoid reflection from 
the boundaries of our box. 



III. APPLICATION: MULTI-PHOTON 
IONIZATION OF NE 

As a first application of our approach, we studied the 
short-pulse, multi-photon ionization of Ne. The elec- 
tric field was taken as linearly polarized. We only ac- 
counted for single-electron excitation and ionization lead- 
ing to ls 2 2s 2 2p 5 n£ bound states or ls 2 2s 2 2p 5 k£ contin- 
uum states in the present, proof-of-principle, calculation. 
Extensions to handle more complex excitations and/or 
double ionization are possible and will be discussed in 
the conclusion of the paper. 

After the wavefunction has been propagated, the rel- 
evant information is extracted by standard projection 
techniques. This requires the ground-state wavefunc- 
tion of the Ne atom, ^>q, obtained either by imaginary 
time propagation or exact diagonalization, and the un- 
perturbed states, *i>® L , where the label 7 represents the 
collection of quantum numbers needed to define the state 



of a bound or free electron in the field of the resid- 
ual atomic ion, asymptotically. In practice, these states 
are constructed from a linear combination of products of 
bound excited states of the Ne + ion coupled to a bound 
or continuum function of the additional electron. 

The quantities of interest for the present work are the 
total survival probability of the initial state, 

P = | < *o|* > | 2 , (19) 
the probability of finding a given (7, L) state, 

P 7 ,lH <*°J*> | 2 , (20) 
and the total probability into a specific L, 

Pl = J2 P ^- (21) 

7 

In computing P 7i l, the time-propagated wavefunction 
is projected onto the singly excited states with an en- 
ergy below (excitation) or above (ionization) the single- 
ionization threshold leading to Ne + . In practice, we com- 
pute the bound-state fraction and get the contribution 
from the continuum by subtraction plus the loss in the 
norm due to the masking function. 

Figure 1 shows our results for the response of a neon 
atom in its ground state (2p 6 ) 1 S to the effect of a 10- 
cycle laser pulse with a sin 2 envelope. Since the frequency 
of the laser pulse is large (sufficient to ionize the atom 
by absorption of a single photon), only a few values for 
the total angular momentum L of the system have to 
be coupled to obtain converged results for single ioniza- 
tion, with relatively small dimensions m of the Krylov 
space. As seen from the figure, the ionization process 
with the highest probability indeed is ionization to the 
P-continuum, i.e., effectively a one-photon absorption 
process leading to a free electron with £ — or £ — 2, re- 
spectively. The survival probability for the ground state 
is approximately 60%, while the probability for excita- 
tion is just under 10%. Finally, the probability for the 
ejected electron to have an orbital angular momentum 
of £ — 1 or I = 3, i.e., forming an L = or L = 2 state 
of the e — Ne + scattering problem, is small but not zero. 

Figure 2 shows the response of the Ne atom to pulses 
with approximately one-half and one-third of the laser 
frequency used in figure 1. In this case, at least two or 
three photons, respectively, need to be absorbed in order 
to ionize the system. A significantly larger number of 
symmetries (we used L max — 6) must be coupled to get 
converged results for these cases. Note that excitation 
rather than ionization appears as the dominating reac- 
tion process for oj = 0.27 a. u. and the laser parameters 
chosen here. 

The dependence of the ionization probability on the 
laser intensity is plotted in figure 3 for laser frequencies 
of 0.425 a.u. and 0.270 a.u., respectively. For peak in- 
tensities in the range 10 13 — 5 x 10 14 W/cm 2 , the slopes 
in the log-log plot (increases of about two or three orders 
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FIG. 1: (Color online). Laser pulse (top panel), ground- 
state survival and excitation probabilities (center) , and single- 
ionization probabilities (bottom) for the interaction of a short 
laser pulse of frequency 0.845 a.u. and peak intensity of 
3.5 x 10 14 W/cm 2 with a neon atom in its (2p°) 1 S ground 
state. Using different sizes m of the Krylov space and num- 
bers of coupled L- values (L max ), we demonstrate in the center 
panel the numerical convergence of our results for the survival 
probability of the ground state as a function of time. In the 
bottom panel, we show the contributions to the total ioniza- 
tion probability from ionization continua with different values 
of the total orbital angular momentum L. 



of magnitude in the probability per one order of mag- 
nitude increase in the intensity) are consistent with the 
expectation from lowest-order perturbation theory that 
ionization is effectively caused by two-photon or three- 
photon processes, without hitting any resonances. For 
higher laser intensities, the curve flattens because of both 
saturation and double ionization. The description of the 
latter processes is, in principle, also possible with the 
current method. However, it requires the inclusion of 
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FIG. 2: (Color online). Ground-state survival (left scale) and 
total excitation and ionization probabilities (right scale) for 
laser frequencies of 0.425 a.u. and 0.270 a.u. and a peak in- 
tensity of 3.5 x 10 14 W/cm 2 . 
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FIG. 3: (Color online). Ionization probability vs. laser peak 
intensity at laser frequencies of 0.270 a.u. and 0.425 a.u. 



double-continuum states with a Ne 2+ core in the cur- 
rent expansion and, therefore, significantly more compu- 
tational resources. 

As a further check of our present work, we now con- 
sider the generalized cross section for two-photon ioniza- 
tion. Having obtained the total two-photon ionization 
rate by propagating the wavefunction in a longer 
pulse (30 optical cycles in the present case), the gen- 
eralized cross section for two-photon ionization is 
obtained as described by Charalambidis et al. [Uj]. In 
figure 4, we compare our non-perturbative results for a 
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FIG. 4: (Color online). Generalized cross section for two- 
photon ionization of Ne(2p 6 ) 1 5 as a function of photon energy. 
The current results (circles) are compared with the i?-matrix 
Floquet predictions of McKenna and van der Hart [26( ] . 
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FIG. 5: (Color online). Excitation probability of the 
Ne(2p 5 3s) 1 P state for a 30-cycle laser pulse with photon en- 
ergies of 15 eV (right scale), 16 eV (left scale), and 17 eV 
(right scale) and a peak intensity of 2.0 x 10 13 W/cm 2 . 



few photon energies to the i?-matrix Floquet predictions 
of McKenna and van der Hart 26]. We note satisfac- 
tory agreement for energies away from the first resonance 
structure corresponding to the intermediate (2p 5 3s) 1 P° 
state. Since the Floquet approach effectively corresponds 
to an infinitely long "pulse" and hence a sharp photon en- 
ergy, it can resolve this structure, while we get a broader 
maximum due to the frequency width of our pulse. The 
shift in the energy position of the resonance is due to the 
different structure models used in the two calculations. 

To further illustrate the effect of the resonant 
(2p 5 3s) 1 P state (around 16 eV in our model) we present 
the excitation probability for three photon energies in 
figure 5. Away from the resonance, at 15 eV and 17 eV, 
there are several Rabi oscillations between the ground 
state and the excited state during a 30-cycle pulse, and 
the maximum probability for excitation is about 5% dur- 
ing these oscillations (right scale of figure 5). On the 
other hand, we just reach the first maximum in the exci- 
tation probability for the resonance energy of 16 eV after 
30 cycles, and the value of that maximum is above 90% 
(left scale of figure 5). This shows the strong effect of the 



energy detuning on the frequency and the amplitude of 
the Rabi oscillations. 



IV. CONCLUSIONS AND OUTLOOK 

We have described a general method for treating the 
interaction of a strong attosecond laser pulse with a com- 
plex atom. The approach combines a highly flexible 
B-spline i?-matrix method for the description of the ini- 
tial state, other bound states in the system, the ionic 
core, and the interaction of the free electron with the 
residual ion after ionization, with an efficient Arnoldi- 
Lanczos scheme for the time propagation of the TDSE. 
The major advantages of the method are i) its generality 
and applicability to any complex many-electron target, 
ii) the possibility of generating highly accurate target 
and continuum descriptions with relatively small config- 
uration interaction expansions, and iii) an efficient time- 
propagation technique. In the current paper we limited 
the continuum states to include only singly ionized states. 
To extend this to doubly ionized targets requires that 
we allow two i?-matrix orbitals outside a doubly charged 
ionic core. In principle, this is a straightforward exten- 
sion of the current codes, but in practice the size of the 
matrix blocks will increase dramatically. The critical ad- 
vantage of the present approach is that non-orthogonal 
basis sets should significantly reduce the size of the con- 
figuration expansion compared to approaches based on 
orthogonal sets. 

In the future, we plan to further analyze and improve 
the numerical efficiency of the method, particularly by in- 
vestigating different schemes of setting up the matrices. 
Prime candidates are expansions in other complete bases 
such as finite-element discrete-variable representations. 
The use of many-electron expansions in non-orthogonal 
basis sets also necessitates developing efficient, new ap- 
proaches to the time propagation of the wavefunction. 
While we have described one possibility in this paper, it 
is not the only one and likely far from the best approach. 
We are actively investigating other methods, using ap- 
proximate and easily computed inverses, which do not 
require any factorization or diagonalization of large ma- 
trices. Most critically, we need to extend the current BSR 
method to treat two free electrons outside an ionic core, if 
we are to treat problems involving multi-photon double- 
ionization of complex targets and compare with recent 
free-electron laser experiments such as those reported in 
refs. [13, Hl| • Finally, looking at angle-differential observ- 
ables will also require a reliable method to extract the 
relevant information, such as amplitudes from single and 
double ionization, from the time-propagated wavefunc- 
tion . All of these issues are currently under active 
investigation by our collaboration. 
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